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ABSTRACT 

Rapid orbital drift of macroscopic dust particles is one of the major obstacles against planetesimal formation 
in protoplanetary disks. We reexamine this problem by considering porosity evolution of dust aggregates. We 
apply a porosity model based on recent A^-body simulations of aggregate collisions, which allows us to study 
the porosity change upon collision for a wide range of impact energies. As a first step, we neglect collisional 
fragmentation and instead focus on dust evolution outside the snow line, where the fragmentation has been 
suggested to be less significant than inside the snow line because of a high sticking efficiency of icy particles. 
We show that dust particles can evolve into highly porous aggregates (with internal densities of much less 
than 0.1 gem -3 ) even if collisional compression is taken into account. We also show that the high porosity 
triggers significant acceleration in collisional growth. This acceleration is a natural consequence of particles' 
aerodynamical property at low Knudsen numbers, i.e., at particle radii larger than the mean free path of the gas 
molecules. Thanks to this rapid growth, the highly porous aggregates are found to overcome the radial drift 
barrier at orbital radii less than 10 AU (assuming the minimum-mass solar nebula model). This suggests that, if 
collisional fragmentation is truly insignificant, formation of icy planetesimals is possible via direct collisional 
growth of submicron-sized icy particles. 

Subject headings: dust, extinction — planets and satellites: formation — protoplanetary disks 



1. INTRODUCTION 

Growth of dust particles is a key process in protoplanetary 
disks. Current theories of planet formation assume kilometer- 
sized solid bodies called "planetesimals" to form from dust 
contained in protoplanetary disks. Being the dominant com- 
ponent of disk opacity, dust also affects the temperature and 
observational appearance of the disks. Furthermore, dust par- 
ticles are known to efficiently capture ionized gas particles in 
the gas disk, thereby controllin g magnetohydrodynamical be- 
haviors of it (ISano et al.ll2000l) . 

Theoretically, however, it is poorly understood how the dust 
particles evolve into planetesimals. One of the most seri- 
ous obstacles is the radial inward drift of macroscopic aggre- 
gates due to the gas dr ag (IWhipplelll972llAdachi et al.ll 19761: 
IWeidenschillin gill 9771) . Because of the gas pressure support 
in addition to the centrifugal force, protoplanetary disks tend 
to rotate at sub-Keplerian velocities. By contrast, dust parti- 
cles are free from the pressure support, and hence they tend 
to rotate faster than the gas disk. The resulting head wind act- 
ing on the dust particles extracts their angular momentum and 
thus causes their drift motion toward the central star. In or- 
der to go beyond this "radial drift barrier," dust particles must 
decouple from the gas drag (i.e., grow l arge) faster than they 
drift inward. However, previous work by Brauer e t al.l(t2008al) 
showed that dust particles finally fall onto the central star un- 
less the initial dust-to-gas mass ratio is considerably higher 
than the canonical interstellar value. 

Several mechanisms have been raised so far regarding how 
dust particles overcome the radial drift barrier. A clas- 
sical idea is that dust particles "jump" across the barrier 
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by forming a gravitationally unstable thin dust layer at the 
midp lane and directly collapsing into planetesimal-size ob- 
jects dSafronovfi969t iGoldreich & Ward|[l97l lHavashi et all 
QUI. However, this classical scenario has been challenged 
by the fa ct that the dust layers are easily stirred up by disk tur - 
bulence (IWeidenschilling & Cuzzilll993t iTurner et all 120101) . 
Moreover, the dust sublayer is known to induce the Kelvin- 
Helmholtz instability, which prevents further sedimentation 
of dust even without disk turbulence unless the dus t -to-ga s 
surface density ratio is considerably high (ISekiyal 11998b . 
Recently, a two-fluid instability of du st and gas has been 
discovered (lYoudin & Goodman! 120051) . which can lead to 
fast formation of gravitationally bound dust clumps (e.g., 
Johansen & Youdinll2007l Uohansen et al J 120071: iBai & Stand 



2010al) . However, this mechanism requires marginally de- 



coupled dust particles, the formation of which is already 
questioned by the radial drift barrier itself. Other possibil- 
ities include the tra pping of dust particles in vortices (e.g., 
iBarge & Sommerial Tl99l [Klahr & Hennind 119971) and at 
gas pressure maxima (e.g . . iKretke & Linl 120071: iBrauer et al.l 
120081* [Suzuki et al.1120 1 Ob IPinilla et al.ll2012l) . 

This study reexamines this problem by considering a new 
physical effect: porosity evolution of dust aggregates . Mos t 
previous coagulation models (e.g.. iNakagawa et al.l 11981 



iTanaka et aT1l2005b IBrauer et al.ll2008at iBirnstiel et al.ll20T 
assumed that dust particles grow with a fixed internal den- 
sity. In reality, however, the internal density of aggregates 
does change upon collision depending on the impact energy. 
The evolution of porosity directly affects the growth history of 
the aggregates since the porosity det ermines the couplin g of 
them to the gas m otion. For example. lOrmel et al.l (120071) and 
IZsom et ail d201 ll) simulated dust growth with porosity evolu- 
tion at fixed disk orbital radii and found that porous evolution 
delays the settling of dust onto the disk midplane. However, 
how the porosity evolution affects the radial drift barrier has 
been unaddressed so far. 
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It has been studied over last two decades how the 
intern al structure changes upon collision by laboratory 
(e.g., iBlum & Wurml |2000t IWeidling et al.l 120091) and nu- 



merical (e.g. IDominik & Tielensl 119971: IWada et afl l200l 
iSuvamaet al.l 1 120081 120121) collision experiments. One ro- 
bust finding of these studies is that aggregates grow into 
low-density, fractal objects if the impact energy is lower 
than a threshold E m \ i determined by ma t erial properties 
dBfum & Wurml 120001: ISuvama et all 120081: lOkuzumi et all 
120091) . The fractal dimension df of the resulting aggregates 
depends weakly on the size ratio between targets and pro- 
jectiles, and falls below two when the target and projectil e 
have similar sizes dMukai et al.lll992t lOkuzumi et al.ll2009l) . 
The fractal dimension of two is equivalent to an internal den- 
sity decreasing inversely proportional to the aggregate radius. 
The density decrease occurs because each merger event in- 
volves the creation of "voids" whose v olume is comparable to 
those of the aggregat es before merger (lOkuzumi et al.ll2009l) . 
ISuvama et alj (120081) estimated the collision energy of ag- 
gregates in protoplanetary disks as a function of size, and 
showed that aggregates composed of 0.1/mi sized particles 
undergo fractal growth in planet-forming regions until their 
size reaches centimeters. This means that the building blocks 
of planetesimals should have once evolved into very fluffy ob- 
jects with mean internal densities many orders of magnitude 
lower than the solid material density. 

More strikingly, recent A^-body experiments suggest that 
the porosity of aggregates can be kept considerably high 
even after t he col lision energy exceeds the threshold E m \\. 
IWada et af] (120081) numerically simulated head-on collisions 
between equal-sized fractal aggregates of df « 2 and found 
that the fractal dimension after the collision does not ex ceed 
2.5 even at high collision energies. ISuvama et all d2008b con- 
firmed this by repeating head-on collisions of the resulting 
aggregates at fixed collision velocities. Furthermore, com- 
paction is even less efficient in offset collisions, where the col- 
lision energy is spent f or stretching rather than compaction of 
the m erged aggregate (IWada et alj|2007t iPaszun & Dominikl 
120091) . These results mean that the creation of voids upon 
merger is nonnegligible even when the impact energy is large; 
in other words, the voids are only imperfectly crushed in colli- 
sional compaction. Because of technical difficulty, these theo- 
retical predictions have not yet been well tested by laboratory 
nor microgravity experiments. Nevertheless, it is worth inves- 
tigating how aggregates grow and drift inward if they evolve 
into highly porous objects. 

In this study, we simulate the temporal evolution of the 
radial size distribution of aggre gates using the adve ction- 
coagulation model developed by iBrauer et al.l (12008 al) . Un- 
like the previous work, we allow the porosities of aggregates 
to change upon collision, depending on their impact energies. 
To do so, we adop t the "v olume-averaging method" proposed 
by lOkuzumi et al.l (120091) . In this method, aggregates of equal 
mass are regarded as having the same volume (or equivalently, 
the same internal density) at each orbital distance, and the 
advection-coagulation equation for the averaged volume is 
solved simultaneously with that for the radial size distribution. 
To determine the porosity change upo n collisional st i cking , 
we use an analytic recipe presented by ISuvama et al.l (120121) 
that well rep roduces the collis i on outcomes of rece nt jV-body 
simulations (IWada et al.l2008t ISuvama et al.ll2008l) as a func- 
tion of the impact energy. These theoretical tools allow us to 
study for the first time how the porosity evolution affects the 
growth and radial drift of dust aggregates in protoplanetary 



disks. 

In order to clarify the role of porosity evolution, we ig 
nore many other effects rele vant to aggrega t e collision, in 



eluding Coulomb interac tion (|pkuzumill2009|: lOkuzumi et al.1 
201 la! bt iMatthews et alJl2012h . bouncing (IZsom et al.ll2010L 



201 lb IWindmark et al.ll2012l). and collisional fra gmentation 
(IBrauer et al.ll2008allbt iBirnstiel et al.ll2009l 120121) . Coulomb 
repulsion due to negative charging can significantly slow 
down the initial fractal growth, but may be negligible once 
the col lisional compaction becomes effective (lOkuzumi et al.l 
1201 lbl) . Bouncing is often observed in laboratory experiments 
for relatively compact (filling factor > 0.1) aggregates, but 
is less likely to occur w hen aggregates are highly porous as 
we co nsider in this study (lLangkowski et al.l2 008; Wad a et ail 
1201 ll) . Seemingly more problematic is fragmentation at high- 
speed collisions. This is particularly so when the aggregates 
are mainly composed of silicate particles, for which catas- 
trophic r ^fisnr£tion_beginsj^ low as a few 
m s' 1 (IBlum & Wurm ll2Q08t IWada et al.ll2009l: iGiittler et all 
120101) . By contrast, collisional fragmentation may be less 
problematic for aggregates made of icy particl es, for which a 
higher sticking threshold has been anticipated (O iokshi et al. 
119931: IDominik & Tielensl 1 19971: iGuiidlach et al.1 1201 ll). For 
instance, A^-body collision experiments by lWada et al.l (120091) 
suggest that aggregates made of 0. 1 fim sized icy grains do not 
experience catastrophic disruption at collision velocities up to 
35-70 m s _1 . For this reason, instead of neglecting collisional 
fragmentation, we focus on dust evolution outside the snow 
line in protoplanetary disks. A more comprehensive model 
including the above mentioned effects will be presented in fu- 
ture work. 

We will show that dust particles evolve into highly porous 
aggregates even if collisional compaction is taken into ac- 
count. Furthermore, we will show that the porosity evolution 
triggers significant acceleration in collisional growth at early 
stages, allowing them to grow across the radial drift barrier 
in inner regions of protoplanetary disks. Interestingly, this 
acceleration involves neither enhancement of the collision ve- 
locity nor suppression of the radial drift speed of marginally 
decoupled aggregates. As we will see, this acceleration is a 
natural consequence of particles' aerodynamical property at 
low Knudsen numbers, i.e., at particle radii larger than the 
mean free path of the gas molecules, and the porosity evolu- 
tion only allows the dust aggregates to reach that stage with 
small aggregate masses. Our model calculation shows that the 
breakthrough of the radial drift barrier can occur in "planet- 
forming" regions, i.e., at < 10 AU from the central star. This 
result suggests that, if the fragmentation of icy aggregates is 
truly negligible, the formation of icy planetesimals is possible 
via direct collisional growth of dust particles even without an 
enhancement of the initial dust-to-gas mass ratio. 

This paper is organized as follows. In Section [2] we de- 
scribe the disk and collision models that we use in this study. 
Simulation results are presented in Section [3] which we in- 
terpret in terms of the timescales for collisional growth and 
radial inward drift in Section |4] The validity and limitations 
of our model are discussed in Section [5] and our conclusions 
are presented in Section [6] 

2. MODEL 

2.1. Disk Model 

We adopt the mi nimum-mass solar nebula (MMSN) model 
of iHavashil dl981l) with a solar-mass central star. The ra- 
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dial profiles of the gas surface density S 9 and disk temper- 
ature T are given by Y, g = 152(r/5 AU)~ 3 ^ 2 gem -2 and 
T = 125(r/5 AU) -1 ^ 2 K, respectively, where r is the distance 
from the central star. In this study, we focus on dust evolution 
outside the snow line, which is located at r ~ 3 AU in the 
adopted disk model. The vertical structure is assumed to be 
in hydrostatic equilibrium, and hence the vertical structure of 
the gas density p g is given by p g = (L g / y2nh g ) exp(-z 2 /2/i 2 ), 
where h g = c 5 /Q is the gas scale height, c s is the isother- 
mal sound speed, and f2 is the Kepler frequency. The isother- 
mal sound speed is given by c s = JkftT/m g , where is 
the Boltzmann constant and m g is the mean molecular mass. 
We assume the mean molecular weight of 2.34, which gives 
m g = 3.9 x 1(T 24 g and c s = 6.7 x 10 4 (r/5 AU)~ l/4 cm s _1 . 

The assumed stellar mass (1 M ) leads to Q = y/GM e /P = 
1.8 x KTV/5 AU)- 3/2 rad s" 1 and h g /r = 0.05(r/5 AU) 1/4 , 
where G is the gravitational constant. 

In reality, protoplanetary disks can be h eavier than the 
MMSN. The gravitational stability criterion dToomre 111964 
Y, g < Qc JttG ~ 5.6 x 10 3 (r/5 AU)~ 7/4 g enr 2 allows the sur- 
face density to be up to about 10 times higher than the MMSN 
value. The dependence of our result on the disk mass will be 
analytically discussed in Section|4] 

Initial dust particles are modeled as compact spheres of 
equal size ao = 0.1 fim and equal internal density pa = 
1.4 g cirT 3 , distributed in the disk with a constant dust-to- 
gas surface density ratio zZ d /zZ g = 0.01. The mass of each 
initial particle is ma = {4jt/3)pocIq = 5.9 x 10~ 15 g. In 
the following, we will refer to the initial dust particles as 
"monomers." We define the radius of a porous aggregate as 
a = [(5/6A0Z£i Z/Lita""*/) 2 ] 1 ' 2 . where N is the number of 
the constituent monome rs and xt (k — 1,2 , . . . , AO is the po- 
sition of the monomers dMukai et al.lll992l) . This definition is 
in accordance with previou s jV-body experiments (IWada et al.l 
I2008t ISuyama et atl 12008b lOkuzumi et all 120091) which our 
porosity model is based on (see Section l2.3.1b . 

Disk turbulence affects the collision and sedimentation of 
dust particles. To include these effects, we consider gas tur- 
bulence in which the turnover time and mean-squared random 
velocity of the largest turbulent eddies are given by = 
and 8v~ = ar>c 2 s , respectively, where ao is the dimension- 
less parameter characterizing the strength of the turbulence. 
The assumption for ti is based on theoretical anticipation for 
turbulence in Keplerian disks (iDubrulle & yaldettarolll992t 
iFromang & Papaloizoull2006t iJohansen et al.ll2006l) . The dif- 
fusion coefficient for the gas is given by D g = S&l = 

aoc 2 JD.. If the gas diffusion coefficient is of the same order as 
the turbulent viscosity, an is equivalent to the so-called alpha 
parameter of lShakura & Sunyaevl(ll973l) . However, we do not 
consider the viscous evolution of the gas disk for simplicity. 
We adopt ao = 10" 3 in our numerical simulations. A higher 
value of ao would cause catastrophic collisional fragmenta- 
tion of aggr egates, which is not considered in this study (see 
Sectionl531l. 



2.2. Evolutionary Equations 

We solve the evolution of the radial size di stribution of 
dust ag gregates using the method developed by iBrauer et al.l 
d2008al) . This method assumes the balance between sedi- 
mentation and turbulent diffusion of aggregates in the vertical 
direction. Thus, the vertical number density distribution of 



aggregates is given by a Gaussian (N / y/2Kh d ) exp(-z 2 /2/i^), 
where N(r, m) is the column number density of aggregates 
per unit mass and h d (r, m) is the s cale height of aggreg ates at 
orbital radius r and with mass m (IDubrulle et al.lll995l) . This 
approach is valid if the coagulation timescale is longer than 
the settling/diffusion timescale, which i s true except for v ery 
tiny particles with short collision times (Zsom et al. 201 ll) . 

The evolution of the radial size distribution N(r, m) is given 
by the vertic ally integrated adve ction-coagulation equation, 
which reads (IBrauer et al.ll2008ah 

dN(r,m) 1 f" 

= — I K(r, m', m — m')N(r, m')N(r, m — m')dm' 

at 2 Jo 

—M(r, m) I K(r, m, m')N(r, m')dm' 
Jo 



1 d 



[ru r (r, m)N(r, m)] , 



(1) 



where v r is the radial drift velocity and K is the vertically in- 
tegrated collision rate coefficient given by 



K(r, m\,mq) 



Ccoll 



2nh d ^h d 7 J_ 



f 

•J — o> 



Auexp 



\dz. (2) 



Here, cr co n is the collisional cross section, h d ,\ and h d ,2 are 
the scale heights of the colliding aggregates 1 and 2, Au is 
the collision velocity, and h d ^2 = {h^\ + ^|) -1 ^ 2 - As men- 
tioned in Section Q] we neglect electrostatic and gravitational 
interactions between colliding aggregates and assume perfect 
sticking upon collision. Thus, the collisional cross section is 
simply given by <r co u = n(a\ + ai) 1 , where a\ and ai are the 
radii of the colliding aggregates. The val idity of neglecting 
fragmentation will be discussed in Section [53] 

The dust scale height h d in sedimenta tion-diffusion equi- 
librium has been analytically obtained bv lYoudin & Lithwickl 
( 120071) . For turbulence of t L = Q, 1 and D g = a D c 2 /Q, it is 
given by 



D.t s 1 + 2Qt s 
ao 1 + Or s 



-1/2 



(3) 



where t s is the stopping time of the aggregates. We use this 
expression in this study. 
For the stopping time, we use 



/Ep) 



/St) 



3/;/ 



4p g v tb A ' 



a < —A 



mfp> 



4a 



9/1 



f ( t Ep) , a > -A 



(4) 



mfp 



mf p » 



where u± = y/8/7ic s and A m f p are the thermal velocity and 
mean free path of gas particles, respectively, and A is the pro- 
jected area of the aggregate. The mean free path is related to 
the gas density as 

m g 

Amfp = , (5) 

V molPg 

where cr mo i = 2 x 10" 15 cm 2 is the collisional cross sec- 
tion of gas molecules. Our gas disk model gives A m fp = 
120(r/5 AU) 11 / 4 cm at the midplane. Equation (0]i satisfies the 
requirement that the stopping time must obey Epstein's law 



t s 



tf p) at a <s A 



mfp 



and Stokes' law f, 



4 St) at a » A mfp , 



respectively. Since f ( s St) oc af' Ep) , an aggregate growing in the 
Stokes regime decouples from the gas motion more quickly 
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than in the Epstein regime. In reality, Stokes' law breaks 
down when the particle Reynolds number (the Reynolds num- 
ber of flow around the particle) is much greater than unity, but 
we neglect this in our simulation s for simplicity. We will dis- 
cuss this point further in Section I5TT1 
The radial drift velocity is taken as 



2Q,t s 



1 + (Qt s ) 



where 



\v K ) d In r 



(6) 



(7) 



is the ratio of the pressure gradient force to the stellar gravity 
in the radial direction and vk = rQ. is the Kepler velocity 
dAdachi et al.l Il976t iWeidenschillind 1 1 9771 iNakagawa et all 
119861) . The radial drift speed has a maximum tjvk, which is 
realized when Qf s = 1 . In our disk model, 77 scales with r as 
77 = 4.0 x 10~ 3 (r/5 All) 1 ' 2 , and the maximum inward speed 
tjvk = 54 m s _1 is independent of r. Since 77 is proportional 
to the gas temperature, the maximum d rift speed would be 
somewhat lower in co lder disk models dKusaka et al.lll970t 
iHirose & Turner! 1201 lb . Equation (|6j neglects the frictional 
backreaction from dust to gas assuming that the local dust-to- 
gas mass ratio is much lower than unity or the stopping time of 
aggregates dominating the dust mass is much longer th an Q" 1 . 
We examine the validity of this assumption in Section l5.2.1l 

In this paper, we also consider the collisional evolution of 
aggregate porosities. We treat the mean volume V = (4n/3)a 3 
of aggregates with orbital radius r and mass m as a time- 
dependent quantity. The evolutionary equation for V(r, m) is 
given by 



d(VN) 1 



[Vi+ 2 K](r, in' , m — rri) 



dt i. jo 

xN(r, m')N(r, m — m')dm' 

K(r, m, m')N(r, m')dm 



1 d 

— — [rv r (r, tn)V(r, m)N(r, m)] , (8) 
r or 



where 



[Vi+2K](r,m u m 2 ) = 



2nhd.\h 



Vi+2A(;exp — dz 



(9) 

with Vi+2 being t he volume of merged aggregates (described 
in Section 12.3.1b . Equation (O is id entical to the original 
evolutionary equation for V derived by lOkuzumi et al.l (120091 
their Equation (16)) except that we here take the vertical inte- 
gration of the equation and take into account the radial advec- 
tion of dust. In deriving Equation ©, we have assumed that 
the dispersi on of the volume is s ufficiently narrow at every r 
and m (see lOkuzumi et all 120091) . This "volume-averaging" 
approximation allows to follow the porosity evolution of ag- 
gregates without solving higher-order moment equations for 
the volume, and hence with small computational costs. This 
approximation is valid unless the porosity distribution at fixed 
r and m is significa ntly broadened by, e. g., collisional frag- 
mentation cascades ( lOkuzumi et al .1120091) . 

2.3. Dust Model 





1+2 



(a) 



(b) 



(c) 



Figure 1. Schematic illustration of our porosity change model, (a) Porous 
aggregates with volumes V[ and V2 before contact, (b) Just after contact. 
At this moment, the volume of the new aggregate is given by Vi+2,hs = 
Vi + V2 + V vo [ L [, whe re V vo id = ^void(^l> ^2) is the volume of newly formed 
voids (Equation 1111 ). If the collision energy £[ mp is much smaller than the 
rolling energy E ro \i, the final volume of the new aggregate is equal to Vi+2,hs- 
(c) If £j mp > E m n, collisional compression occurs. In this case, the final 
volume Vi+2(< Vi+2,hs) depends on £ imp . 



2.3.1. Porosity Change Recipe 

The functional form of Vi+2 determines the evolution of ag- 
gregate porosities in our simulation. In this study, we give 
V\+2 as a function of the volumes of the colliding aggregates, 
V\ = V(r,mi) and V2 = V(r,m 2 ), and the impact energy 



J imp 



«i 1 ?Ti2 At; 2 / [2(»j 1 + ni2)]. Before introducing the fi- 
nal form of our porosity change recipe (Equation (TT5t ), we 
briefly review recent A^-body collision experiments which our 
recipe is based on. 

Collisional compression depen ds on the ratio between 
Ej m „ and the "rolling e nergy" E m u ( Dominik & Tielenslll997l : 
iBlum & Wurml l2000t IWada et al.l 120071) . The rolling en- 
ergy is defined as the energy needed for one monomer to 
roll over 90° on the surfa ce of another monomer in contact 
(iDominik & Tielenslll997l) . When Ei mp <sc isVoii, two aggre- 
gates stick without visible restructuring (so-called hit-and- 
stick collision; see Figure [T{b)). In this case, the volume of 
the merged aggregate is determined in a geometric manner, 
i.e., independently of Ei mp . When Ei mp > E m n, internal re- 
structuring occurs through inelastic ro lling among constituent 
monomers (IDominik & Tielensl[l997l: see also Figure Q2c)). 
In this case, the final volume V\+2 depends on Eimp as well as 
on V\ and V 2 . 

For hit-and -s tick collisions (Ei mp /E m \\ — > 0), 
lOkuzumi et al.l (120091) obtained an empirical formula 

for Vi +2 , 

V 1+2 = Vl + 2,HS = V l + V 2 + Koid, (10) 

where Vi and V2(< Vi) are the volumes of the two colliding 
aggregates, and 



y void = min -10.99- 1.03 lnl 



V1/V2 + I 



6.94 VV; 



(11) 



is the volume of the voids created in the collision (see Fig- 
ure [Tib))- For Vi ~ V 2 , the void volume is approximately 
equal to Vi, and hence the volume of the new aggregate is 
approximately given by Vi + 2 ~ 3Vi. This is equivalent to a 
fractal relation V oc m 3 ^ df , where df w 2 (see Section 4.2.1 of 
lOkuzumi etaill2009l) . 
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In the limit of Ei mp » E ro u and for head-on c o llision 
of equal-sized aggregates {V\ = V2), ISuvama et al.l (12008b 
showed that Vi+2 obeys the relation 



j imp 



r-Vn-2 
J2 6 /5Vi 



P{V)dV. 



(12) 



Here, P = -dEi mv /dV is the dyn amic compression strength 
of the merged aggregate given by (Wada et al. 2008) 



1+2' 



(13) 



where b = 0.15 is a dimensionless fitting parameter, Vq = 
m o/Po = (4?r/3)aQ is the monomer volume, N\+2 = 2m\ /mo is 
the number of monomers contained in the merged aggregate, 
and pi nt = 2/«i/y is the internal density of the merged aggre- 
gate. If we substitute Equation ( TT3l into Equation ( 1121 . we 
obtain the equation that explicitly gives V\ + 2 as a function of 
E mU and Vu 



V 



1+2 



(3/5) 5 £ n 



Nl 2 bE ron V 



10/3 



-3/10 



(14) 



This equation basically expresses the energy balance in colli- 
sional compression, but some caution is needed in interpreting 
it. First, the initial state for the compression is chosen to be 

V = 2 6 ^ 5 Vi, although the volume just after contact is V — 3V\ 
(see above). This is based on the fact that compaction from 

V = 3V, to V = 2 6/5 V, occurs through partia l compression 
of the new voids, which requires little energy (ISuvama et al.l 
120081) . Second, the dynamic compression strength P depends 
on mass iVi+2 as well as on internal density p mt , meaning that 
P is not an intensive variable. This is due to the fact that dy- 
namically compressed parts in the merged aggre gate have a 
fracta l structure with the fractal dimension of 2.5 (Wada et al. 
120081) . In fact, Equations (fT2l - (fT4l i are more naturally de- 
scribed in terms of variables in the 2.5-dimensional space, 
V f oc Pf e x AWW, and P f = -dE imp /dV f (see 
IWada et al.l 12008b ISuvama et all 120081) . An important point 
here is that aggregates become stronger and stronger against 
dynamic compression as they grow because of the N^Q factor 
in P. 

Equations (TTOb and (fl4l express how the volume of the 
merged aggregate is determined in the limits of E[ mp «: E I0 n 
and Ei mp » E I0 u, respectively. To properly take into account 
intermediate cases (E\ mp ~ E m n ), we adopt an upda ted an- 
alytic formula given recently by ISuvama et al.l (120121) . This 
reads 



Vl + 2 = 



J imp 



3bE m u 

(if vlH 



V 



5/6 

1+2.HS 



J imp 



3bE 



roll 



5/6 



HS > Vf /6 + y 2 5/6 and E imp 



vf) 

< 3bE mU ), 



6/5 



(3/5) 5 (£ imp - 3bE roll ) 
N 5 1+2 bE mU V W/3 



{v\ 16 + y 2 5 



2 5/ T 



-3/10 



( if y f + 2,HS > V l 16 + V 2 6 and £ ™ P > 3^£roll), 



(3/5) £"imp 



+ v; 



-10/3 



N S hp 1/10/3 1 '1+2.HS 

(if V 

u± V 1+2,HS 



J roll ' 

,5/6 



-3/10 



< v' /6 + V 5 ' 6 ), 



where /Vi+2 is now defined as (m\ + mi)lm.Q. Note that this 
equation reduces to Equations ( TTOb when E\ mp «: E m \\, and to 
Equati on (fl4l when V\ = V2 and £i mp » £ ro u- ISuvama et alj 
(120121) derived Equation (TT~5T > by taking into account small en- 
ergy loss in the partial compression of the new voids. In ad- 
dition, unlike Equation (fl4l . Equation ( Tl3T > takes into account 
the cases where colliding aggrega tes have different vol umes 
and masses {V\ + V%, m\ + mi). ISuvama et al.l (120121) con- 
firmed that Equation (TT3T > reproduces the results of numerical 
collision experiments within an error of 20% as long as the 
mass ratio OT2/mi(< 1) between the colliding aggregates is 
larger than 1/16. 

We comment on three important caveats regarding our 
porosity recipe. First, Equation (TTSt is still untested for 
cases where colliding aggregates have very different sizes 
(ni2/mi < 1/16). Therefore, the validity of using Equa- 
tion (TT3T > is at present only guaranteed for the case where 
"similar-sized" (m.2jni\ > 1/16) collisions dominate the 
growth of aggr egates. We will carefully check this validity 
in Section 13.21 Second, Equation ( fTBT ) ignores offset colli- 
sions, in which a considerable fraction of th e impact energy is 
spent for stretching rather than compaction (IWada et al.l2007t 
iPaszun & Dominikll2009l) . For this reason, Equation (fT5l l un- 
derestimates the porosity increase upon collision. Third, we 
do not consider non-collisional compression (e.g., static com- 
pression due to gas drag forces), which could contribute to 
the compaction of very large aggregates. We wi ll dis cuss the 
second and third points in more detail in Section IB~4l 

In addition to V, we need to know the projected area A of 
aggregates to calculate the stopping time t s . Unfortunately, a 
naive relation A = na 2 breaks down when the fractal dimen- 
sion of the aggregate is less than 2, since na 1 increases faster 
than mass for this case while A does not. A projected area 
growing faster than mass means a coupling to the gas becom- 
ing stronger and stronger as the aggregate grows, which is 
cle arly unrealistic. To avo id this, we use an empirical formula 
bv lOkuzumi et al.l d20 09) that well reproduces the mean pro- 
jected area A of aggregates with monomer number /V = m/mo 
and radius a for both fractal and compact aggregates. With 
this formula, all aggregates in our simulations are guaranteed 
to decouple from the gas as they grow. We remark that this 
treatment is only relevant to fractal aggregates with df < 2; 
for more compact aggregates, the empirical formula reduces 
to the usual relation A as na 2 . 

The rolling energy £,011 has not been measured so far for 
submicron-sized icy particles, but can be estimated in the fol- 
lowi ng way. It is anticipated by microscopic friction the- 
ory (iDominik & Tielens| [l995) that the critical rolling force 
Froii = E m n/(7rao/2) is a material constant (i.e., .EVoii is propor- 
tional to the monomer radius «o)- Recently, a rolling force of 
F m \\ = (1 .15 ±0.24 ) x 10~ 3 dyn has been m easured for micron- 
sized ice particles (Gu ndlach et al.ll201 ll) . Given that F [0 u is 
independent of ao, the measured force implies the rolling en- 
ergy of E roa = (nao/2)F ma 1.8 x 10~ 8 erg for a = O.lyum. 
We use this value in our simulations. 

2.3.2. Collision Velocity 

We consider Brownian motion, radial and azimuthal drift, 
vertical settling, and turbulence as sources of the collision 
velocity, and give the collision velocity Au as the root sum 
square of these contributions, 

(15) Av = yJ(Av B ) 2 + (At,,.) 2 + (At, ) 2 + (Ay,) 2 + (A<;,) 2 , (16) 
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where Aug, Av r , Ay^, Ay,, and Ay, are the relative velocities in- 
duced by Brownian motion, radial drift, azimuthal drift, ver- 
tical settling, and turbulence, respectively. 
The Brownian-motion-induced velocity is given by 



Av E 



Tlifl \ Wl2 



8(mi + ni2)kftT 



(17) 



where m\ and mi are the masses of the two colliding aggre- 
gates. 

The relative velocity due to radial drift is given by Ay,. = 
- v r (t s ,2)\, where t Si i and t s ,i are the stopping times of 
the colliding aggregates, and v r is the radial velocity given by 
Equation ©. Similarly, the relative velocity due to differen- 
tial azimuthal motion is given by Ay^ = \v'^(t St {) - 1^(^,2)!, 
where 

"< ~ "TT(Q^ (18) 

is the deviation of the azimuthal velocity from the local 
Kepler velocity (lAdachi et alJ fl976t IWeidenschillingl 1 1 977b 
iNakagawa et al.l[l986l) . Here, we have neglected the back- 
reaction from dust to gas as we alre ady did for the radial ve- 
locity (see Sections [2.3.2l and l5.2. 11 1. 

For the differential settling velocity, we assume Ai> z = 
k-(f.u) - v z (t s>2 )\, where 



1 + Qr o 



(19) 



Equation ( [19l i reduces to the terminal settling velocity u z = 
-Cl 2 t s z in the strong coupling limit Of s «; 1, and to the 
amplitude of the ve rtical oscillation velocity at Qf s » 1 
dBrauer et al.ll2008al) . 

For the turbulence-driven relative velocity, we use an 
analytic formu l a for Kolmogorov turbulence derived by 
lOrmel & Cuzzil J2007L their Equation (16)). This analytic for- 
mula has three limiti ng expressions (Equations (27)-(29) of 
lOrmel & Cuzzill2007l) : 



Ay, 



t s ,i •« t„ 



(1.4. 



,1.7) XSVg yfUiZl, 

1 1 



1 + Qr v 



1 + £if t 



1/2 



t n «: f,,i «; £2 
Q.t s i » 1 , 



6v g , which is reached 



(20) 

— 1/2 

where Re, is the turbulent Reynolds number, t n = Re, ' ti 
is the turnover time of the smallest eddies, and the numeri- 
cal prefactor (1.4 ... 1.7) in the second equality depends on 
the ratio between the stopping times, t s ,2/t s ,i- The turbu- 
lent Reynolds number is given by Re, = D g /v mo \, where 
Vmoi = tWmfp/2 is the molecular viscosity. For t s ,\ ~ t s> 2, 
the maximum induced velocity is Ay, 
when Of 5 1 as 1 . 

When two colliding aggregates belong to the Epstein 
regime and their stopping times are much shorter than t n (<K 
Q~'), the relative velocity driven by sedimentation and turbu- 
lence is approximately proportional to the difference between 
the mass-to-area ratios m/A of two colliding aggre gates. In 
this case, as pointed out by lOkuzumi et al.l (1201 lal) . the dis- 
persion of the mass-to-area ratio becomes important for frac- 
tal aggregates of dj < 2, since the mean mass-to-area ratio 
of the aggregates approaches to a constant and hence the dif- 
ference in m/A(m) vanishes. To take into account the disper- 
sion effect, we evaluate the differential mass-to-area ratio as 




orbital radius r [AU] 



10' 10 2 
orbital radius r [AU] 



Figure 2. Aggregate size distribution AS^/A log m at different times / for the 
compact aggregation model (pj nt = 1 .4 g cirT 3 ) as a function of orbital radius 
r and aggregate mass m. The dotted lines mark the aggregate size at which 
Clt s exceeds 0.1. 

\A(m/Af = \m x [A x - m 2 /A 2 \ 2 + e^mJAtf + (m 2 /A 2 ) 2 ], 
where Aj = A(rrij) (j = 1,2) are the mean p rojected area 
of aggregates with mass mj (see Section 12.3. U and e is the 
standard deviation of the mass-to-area ratio divided by the 
mean ( for the derivation, see the Appendix of lOkuzumi et al.l 
1201 lah . We assume e — 0. 1 in accordance with the numerical 
estimate bv lOkuzumi et alJ d201 lal) . 

2.4. Numerical Method 

We solve Equations (HJ and ([8]) numerically with an explicit 
time-integration scheme and a fixed-bin method. The radial 
domain is taken to be outside the snow line, 3 AU < r < 
150 AU, discretized into 100 rings with an equal logarithmic 
width Aln(r[AU]) = (In 150-ln3)/100. The advection terms 
are calculated by the spatially first-order upwind scheme. We 
impose the outflow and zero-flux boundary conditions at the 
innermost and outermost radii (r = 3 AU and 150 AU), re- 
spectively; thus, the total dust mass inside the domain is a 
decreasing function of time. Our numerical results are un- 
affected by the choice of the boundary condition at the out- 
ermost radius, since dust growth at this location is too slow 
to cause appreciable radial drift within the calculated time. 
The coagulation terms are calculated by the method given by 
lOkuzumi et al.l (120091) . Specifically, at the center of each ra- 
dial ring we divide the mass coordinate into linearly spaced 
bins ntk = hno (k — 1,2,..., A^) for m < A^/'«o and loga- 
rithmically spaced bins nik = ntk-\ lO 1 ^ (k = Nbd + 1 , ■ • ■) for 
m > Nhdi nyi, where Nhd is an int eger. We adopt Nf,d = 40; as 
shown bv lOkuzumi et alJ (120091) . the calculation results well 
converge as long as Nbd 40. The time increment At is ad- 
justed at every time step so that the fractional decreases in 
N and VN fall below 0.5 (i.e., At < -0.5(3 In N/dt)- 1 and 
At < -0.5(3 In VAf/df)- 1 ) at all bins. 

3. RESULTS 

3.1. Compact Aggregation 

To begin with, we show the result of compact aggregation. 
In this simulation, we fixed the internal density p mt = m/V of 
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Figure 3. Radial profiles of the total dust surface density at different times 
for the compact aggregation model (pj llt = 1.4 g cm -3 ). 
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Figure 4. Aggregate size distribution Al.j /A log m at r = 5 AU and t = 
2000 yr-^1470 yr for the compact aggregation model. The dashed and solid 
arrows indicate the aggregate sizes at which a = /l m f p and flf, = 1 , respec- 
tively. Shown at the top of the panel is the aggregate ra dius a. The vertical 
bars indicate the weighted average mass (m) m (Equation 121 1 ). 

the aggregates to the material density po = 1-4 g cm 4 , and 
solved only the evolutionary equation for the radial size dis- 
tributi on N(r, m) (Equatio n (fl]i), as done in previous studies 
(e.g.. iBrauer et al.ll2008al) . Figure [2] shows the snapshots of 
the radial size distribution at different times. Here, the distri- 
bution is represented by the dust surface density per decade 
of aggregate mass, AE^/Alogm = ln(10)m 2 7V(r, m). At each 
orbital radius, dust growth proceeds without significant radial 
drift until the stopping time of the aggregates reaches Qf s ~ 
0.1 (the dotted lines in Figure|2]i. However, as the aggregates 
grow, the radial drift becomes faster and faster, and further 
growth becomes limited only along the line £lt s ~ 0.1 on the 
r—m plane. Figure [3] shows the evolution of the total dust sur- 
face density = J mNdm = f(AX^/Alogm)dlogm. We 
see that a significant amount of dust has been lost from the 
planet-forming region r < 30 AU within 10 5 yr. In this re- 
gion, the dust surface density scales as r~'Q and hence the 



4 It can be analytically shown iBirnstiel et al. 2012) that the dust surface 
density profile obeys a scaling l.j oc -^Z g /(r 2 fi) (oc r _1 for Y. q oc r~ 3 ' 2 ) when 



dust-to-gas surface density ratio oc r~'/£ 9 K r ^ 2 decreases 
toward the central star. 

Figure [4] shows the evolution of the dust size distribution 
observed at r — 5 AU. Here, in order to characterize the 
typical aggregate size at each evolutionary stage, we introduce 
the weighted average mass {m) m defined by 



(m)„ 



J m 2 Ndm 
J mNdm 



m— d log m. 



A log m 



(21) 



The weighted average mass approximately corresponds to 
the aggrega te mas s at the peak of AX,/ / A log w (see, e.g., 
lOrmel et al.l 12007b lOkuzumi et al.ll201 lal) . In Figure g] the 
weighted average mass at each time is indicated by the short 
vertical line. At r — 5 AU, the growth-drift equilibrium is 
reached at t « 4000 yr, and the typical size of the aggregates 
is (m) m x 500 g in mass (« 4 cm in radius, » 0.07Q -1 in 
stopping time). Note that the final aggregate radius is much 
smaller than the mean free path A m f p of gas molecules (the 
dashed arrow in Figure|4}, which means that the gas drag onto 
the aggregates is determined by Epstein's law. As we will see 
in the following, porosity evolution allows aggregates to reach 
the Stokes drag regime at much smaller f2f s . 

3.2. Porous Aggregation 

Now we show how porosity evolution affects dust evolu- 
tion. Here, we solve the evolutionary equation for V(r,m) 
(Equation ([8j) simultaneously with that for N(r,m) (Equa- 
tion ([TJi). The result is shown in Figure [5] which displays the 
snapshots of the aggregate size distribution AZ^/Alogm and 
internal density p; nt = m/V at different times f as a function 
of r and m. The evolution of the total dust surface density X</ 
is shown in Figure [6] 

The left four panels of Figure [5] show how the radial size 
distribution evolves in the porous aggregation. At t < 10 3 yr, 
the evolution is qualitatively similar to that in compact aggre- 
gation (Section lTTt . However, in later stages, the evolution is 
significantly different. We observe that aggregates in the in- 
ner region of the disk (r < 10 AU) undergo rapid growth and 
eventually overcome the radial drift barrier lying at £lf , ~ 1 
(dashed lines in Figure [3j within t ~ 10 4 yr. At this stage, 
the radial profile of the total dust surface is hardly changed 
from the initial profile, as is seen in Figure |6] In the outer 
region (r > 10 AU), aggregates do drift inward before they 
reach Q.t , ~ 1 as in the compact aggregation model. However, 
unlike in the compact aggregation, the inward drift results in 
the pileup of dust materials in the inner region (r « 4-9 AU) 
rather than the loss of them from outside the snow line (see 
Figure |6). This occurs because most of the drifting aggre- 
gates get captured by aggregates that have already overcome 
the drift barrier. As a result of this, the dust-to-gas mass ra- 
tio in the inner regions is enhanced by a factor of several in 
10 5 yr. 

The right four panels of Figure |5] show the evolution of the 
internal density p mt = m/V as a function of r and m. First 
thing to note is that the dust particles grow into low-density 
objects at every location until their internal density reaches 
Pint ~ 10 5 — 10~ 3 g cm 4 . In this stage, the internal density 
decreases as p mt ss (m/mo) _1 ' 2 po> meaning that the dust parti- 
cles grow into fractal aggregat es with the fractal dimension 
df » 2 (lOkuzumi et al.l 120091) . The fractal growth gener- 



ally occurs in early growth stages where the impact energy 

radial drift balances with turbulence-driven growth. 
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Porous Aggregation Model 
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Figure 5. Aggregate size distribution AX^/Alogm (left four panels) and internal density p mt = m/V (right four panels) at different times 
aggregation model as a function of orbital radius r and aggregate mass m. The dashed lines mark the aggregate size at which Q.t s exceeds unity. 
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Figure 6. Radial profiles of the total dust surface density T.j at different times 
for the porous aggregation model. 



is too low to caus e collisional compression, i.e., E mri <si E m \\ 
(e.g.. lBlumll20()4t lOrmel et1iI1l2007t TZsom et al.ll2010h . At 
m ~ 10~ 4 — 10 6 g, the fractal growth stage terminates, fol- 
lowed by the stage where collisional compression becomes 
effective (E[ mp » isVoii)- In this late stage, the internal density 
decreases more slowly or is kept at a constant value depending 
on the orbital radius. We will examine the density evolution 
in more detail in Section l3.2.2l 

Figure Q shows the evolution of the mass distribution func- 
tion at r — 5 AU during t « 1200-2500 yr. The evolution of 
the weighted average mass (m) m is shown in Figure [8] It is 
seen that the acceleration of the growth begins when the ag- 
gregate size a exceeds the mean free path of gas molecules, 
/Imfp (the dashed arrow in Figure [7). This suggests that the 
acceleration is due to the change in the aerodynamical prop- 
erty of the aggregates. At a » /t m f p , the gas drag onto the 
aggregates begins to obey Stokes' law. In the Stokes regime, 
the stopping time t s of aggregates quickly increases with size 
(see Section l2l2V This causes the quick growth of the aggre- 




aggregate mass m [g] 

Figure 7. Aggregate mass distribution AS^/Alogm at r = 5 AU and t = 
1289 yr-2450 yr for the porous aggregation model. The dashed and solid 
arrows indicate the sizes at which a = /t m f p and fir s = 1 , respectively. Shown 
at the top of the panel is the aggregate radius a measured at t = 2 450 yr. The 
vertical bars indicate the weighted average mass (m) m (Equation )2 U ). 



gates since the relative velocity between aggregates increases 
with t s (as long as Qt s < 1). As a result of the growth accel- 
eration, the aggregates grow from a as A m f p to Of , « 1 within 
300 yr, which is short enough for them to break through the 
radial drift barrier. 

The decrease in the internal density plays an important role 
on the growth acceleration. More precisely, the low inter- 
nal density allows the aggregates to reach a as A m f p at early 
growth stages, i.e., at small Qt s . In fact, the growth acceler- 
ation was not observed in the compact aggregation, since the 
aggregate size is smaller than the mean free path at all Q.t s < 1 
(see Figure HJ. A more rigorous explanation for this will be 
given in Section|4] 

3.2. 1 . Projectile Mass Distribution 
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time / [yr] 

Figure 8. Weighted average mass (m} m (Equation 12 U ) at r = 5 AU as a 
function of time t for the porous aggregation model. Shown at the right of 
the panel is the corresponding aggregate radius a((m) m ). The dashed and 
solid arrows indicate the sizes at which a((m) m ) = i m f p and flf ( ((m) m ) = 1, 
respectively. 



range 0.1m, < m p < m,. In fact, the projectile mass dis- 
tribution integrated over 0.1m, < nip < m, exceeds 50% 
for all the cases presented in Figure [9] This means that the 
growth of aggregates is indeed dominated by collisions with 
similar-sized ones as required by the limitation of our poros- 
ity model. This is basically the consequence of the fact that 
the aggregate mass distribution AE c i / A log m is peaked around 
the target mass m ~ (m),„ (see Figure |7). The mass ratio 
m p /m, at the peak of m p C mt (m p ) reflects the size dependence 
of the turbulence-driven relative velocity Av,, which is the 
main source of the collision velocity in our simulation. At 
t < 2000 yr (<m),„ < 10 3 g), the dominant projectile mass is 
lower than m,(= (m) m ), since both the target and projectiles 
are tightly coupled to turbulence (i.e., t s {m,), t s (m p ) < t n ) and 
hence Ay, vanishes at equal-sized collisions (see the first ex- 
pression of Equation (|20V). At t > 2000 yr «m)„, > 10 3 yr), 
the target decouples from small turbulent eddies (f. s (m,) > t n ). 
This results in the shift of the dominant collision mode to 
trip ~ m, because At;, no more vanishes at equal-sized col- 
lisions (see the second line of Equation d20ii). 
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Figure 9. Normalized projectile mass distribution per unit logarithmic pro- 
jectile mass, m p C m ,(mp), for a target with mass m, = (m),„ at different times 
/(= 1289 yr-2450 yr) at r = 5 AU for the porous aggregation model (see 
Equation 1221 for the definition of C m , (m p )). The filled circles show the val- 
ues for equal-sized collisions, m„ = m,(= (m) m ). The dotted and solid arrows 
indicate the target mass at which f , = f, ; and Clt s = 1 , respectively. 



As noted in Section 12.3.11 our porosity change model has 
only been tested for collisions between similar-sized aggre- 
gates. To check the validity of using this m odel, we intro- 
duce the projectile mass distribution function (lOkuzumi et al.l 
I2009T) : 



C m ,(m p ) = 



m p K(m p , m,)N(m p ) 
J m' p K(m' p ,m t )N{m' p )dm' p 



m p < m,, (22) 



which is normalized so that J C mi (m p )dm p — 1. The denom- 
inator of C mt (m p ) is equal to th e growth rate t el l o w = dlnmjdt 
of a target having mass m, (see lOkuzumi et al1l2009l) . Hence, 
the quantity C m ,{m p )dm p measures the contribution of projec- 
tiles within mass range [m p , m p + dm p ] to the growth of the 
target. 

Figure [9] shows the projectile mass distribution per unit 
lnm p , m p C mt (m p ), for targets with mass m, = (m) m at r — 
5 AU and at different t. We see that the growth of the 
m, = {m) m target is dominated by projectiles within a mass 



3.2.2. Density Evolution History 

To see the density evolution history in detail, we plot in 
Figure [10] the temporal evolution of the weighted average 
mass {m),„ and the internal density of aggregates with mass 
m = (m),„ at orbital radii r — 5 AU and 20 AU. 

As mentioned above, dust particles initially grow into frac- 
tal aggregates of df ~ 2 until the impact energy E\ mp becomes 
comparable to the rolling energy E m \\. With this fact, one 
can analytically estimate the aggregate size at which the frac- 
tal growth terminates. Our simulation shows that the colli- 
sion velocity between the fractal aggregates is approximately 
given by the turbulence-driven velocity in the strong-coupling 
limit (Equation ( 1201 with t s <k t n ). Assuming that the colli- 
sions mainly occur between aggregates of similar sizes (see 
Section [3.2.1l ). the reduced mass and the collision velocity are 
roughly given by m/2 and SVgRs 1 / Clt s , respectively. In ad- 
dition, we use the fact that fractal aggregates with df * 2 
have the mass-to-area ratio comparable to their constituent 
monomers. This means that the stopping time of the aggre- 
gates is as short as the monomers and is hence given by Ep- 
stein's law. Thus, the impact energy is approximated as 



m A 2 3 / SVgRe 

E\mn ~ — AlC — m 

p 4 ' 8 



1/4, 



fM /3m 
PgVth j \4A 



(23) 



Furthermore, using the definitions for p g , and Re,, we 
have p g v th = (2/^)Z 9 Q and Re, = a D T.gO- mol /(2m g ) for the 
midplane. Substituting them into Equation ( 1231 and using 



dv g = tJodc s and m/A » mol(nafy = 4po«o/3, we obtain 



J imp 



3tt 2 
32 V2 



3/2 2 



1/2 



(24) 



Thus, the impact energy is proportional to the mass. We define 



the rolling mass m m \\ by the condition E{ mv 
Equation d24b , the rolling mass is evaluated as 



Iroll- 



Using 



mrou ' 



32 V2 grou / 



3tt 2 cia^yZgO-moi) VPoflo 



1/2 



I. 



10- 



lO- 3 



-3/2 



{— 

\100K/ \ 100 gem 



,3/2 



i-2 
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Figure 10. Temporal evolution of the weighted average mass (m) m and the 
internal density Pi nl ((m),„) at orbital radii r = 5 AU (upper panel) and 20 AU 
(lower panel). Shown at the top of the panels is the aggregate radius a((m) m ) 
at each orbital radius. The triangles, circles, diamonds, and square mark the 
sizes at which E im p = E m \\, a = /t m f p , t s = t r] , and Slt s = 1, respectively. At 
r = 20 AU, dust growth stalls due to the radial drift barrier (cross symbol) 
before reaching £lt s = 1 . 



10 3 dyn /\ 1 g cm 



Po 



-2 



do 



0.1 pm 



(25) 



where 
tion|23 



we have used that E m n = (nao/2)F m n (see Sec- 
Using the relations a « (m/mo) ^ 2 Oq and p m « 
(to/too)~^ 2 Po ior df » 2 aggregates, the corresponding radius 
and internal density are found to be 



a roU ~ 1 cm 



'Wroll 

io- 4 g 



1/2 



Pint,roll ~ 10 5 g Cm 3 



»?roll 
10- 4 £ 



-1/2 



(27) 



The triangles in Figure [10] mark the rolling mass at r — 5 AU 
and 20 AU predicted by Equation (|25l l. The analytic predic- 
tion well explains when the decrease in pint terminates. 

The density evolution is more complicated at m > m m \\, 
where collisional compression is no longer negligible (i.e., 
E\ mv > E m \\). At r = 5 AU, the internal density is approx- 



imately constant until the stopping time reaches Qt s = 1, and 
then decreases as p; nt oc m~ l/5 . At r = 20 AU, by contrast, the 
density is kept nearly constant until m ~ 10 2 g (a ~ 10 2 cm), 
and then decreases as p; nt a m 4 ' 8 . 

As shown below, the density histories mentioned above 
can be directly derived from the porosity change recipe we 
adopted. Let us assume again that aggregates grow mainly 
through collisions with similar-sized ones (ni\ as TO2 and 
V\ ~ V2). In this case, the evolution of p mt at E\ mp » E m \\ 
is approximately given by Equation ( fl4l i. Furthermore, we 
neglect the term (2 V^ 6 )~ 4 in Equation ( fl4l ) assuming that the 
impact energy is sufficiently large (which is true as long as 
Qf s < 1; see below). Under these assumptions, the internal 
density of aggregates after collision, p mt = 2toi/Vi+2, is ap- 
proximately given by 



Pint 



3/2 



J imp 



Ni +2 bE mn 



3/10 



(28) 



where N\+2 = 2m! /mo. Since the impact energy E lmp at 
TOi(Ai/) 2 /4 is proportional to A^i + 2(Ay) 2 , Equation d28| i implies 
that 

p in , cc (Ai;) 3/5 TO- 1/5 , (29) 



where we have dropped the subscript for mass for clarity. 
Equation ( 1291 gives the relation between p; nt and to if we 
know how the impact velocity depends on them. Explicitly, 
if At; oc mPp y , Equation d29l leads to 



Pint oc m 



(3/3-l)/(5-3y) 



(30) 



In our simulation, the main source of the relative velocity 
is turbulence. The turbulence-driven velocity depends on t s 
as Av, oc t s at t s «: t n and Au, oc y/T, at t n «: t s <sc t L (- Q~') 
(see Equation (TSUI). As found from Equation ©, the stopping 
time depends on p mt and m as t s oc m/A oc m/a 2 oc m x ^pS in 

the Epstein regime (a <K A m f p ) and as t s oc ma /A oc m 2 ^p l P 
in the Stokes regime (a » /l m f p ). Using these relations with 
Equation d30l ), we find four regimes for density evolution, 



Pint 





l/4 ) 

-1/8^ 





a «: /} mfp and t s 
a » /l m f p and t s 



« tr, 



a <k /l m f p and t n «: t s «: f^, 
a » /l m f p and t n «: f s «: r^. 



(31) 



The circles, diamonds, and square in Figure [10] mark the 



size at which a = A m f p (i.e., t 



(Ep) 



tf\ t, 



t„, and Q? s = 1, 



respectively. At r = 5 AU, the sizes at which a = A m f p and 
f t = t n nearly overlap, and hence only two velocity regimes 

ISO 



t s = t ( s Ep) <K tjj and t n <s f s = fj J '' f L are effectively relevant. 
For both cases, Equation OTb predicts fiat density evolution. 
At r = 20 AU, there is a stage in which t s » and a <K 
(26) -^mfp, for which Equation (f3TT > predicts p; nt oc to -1 ' 8 . These 
predictions are in agreement with what we see in Figure [10] 

Equation ( f28l ) does not apply to the density evolution at 
Q?. s > 1, where the collision velocity no more increases 
and hence collisional compression becomes less and less ef- 
ficient as the aggregates grow. However, if we go back to 



Equation (fl4l i and assume that the impact energy E imp is suf 



imp 

ficiently small, we obtain V\+2 ~ 2 6 ^ 5 Vi, or equivalently 
i+2/«Jj^ 2 * Vi/to^ 5 , where toi + 2 = 2toi is the aggregate 
mass after a collision. This implies that V/m 6 ^ 5 is kept con- 
stant during the growth, i.e., V oc m 6 ' 5 , and hence we have 



V 
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Pint = tn/V °c m '^ 5 . This is consistent with the density evo- 
lution at Qt s > 1 seen in the upper panel of Figure [TOl 

4. CONDITION FOR BREAKING THROUGH THE 
RADIAL DRIFT BARRIER 

In this section, we explain why porous aggregates over- 
come the radial drift barrier in the inner region of the disk. 
We do this by comparing the timescale of aggregate growth 
and radial drift. We assume that dust aggregates grow mainly 
through c ollision s with similar-sized aggregates. As shown 
in Section [3.2.1l this is a good approximation for the growth 
of aggregates dominating the total mass of the system (i.e., 
aggregates with mass m ~ (m) m ). The growth rate of the ag- 
gregate mass m at the midplane is then given by 



dm 2,/ 
-37 = p</cr coll Ai; = — AAv, 
at V2jrh d 



(32) 



where p<j = Z f //( yfhrhd) is the spatial dust density at the 
midplane, and we have approximated <j co \\ as the projected 
area A. Equation ( 1321 can be rewritten in terms of the growth 
timescale as 



' d In m 

i dt 



hd ml A 4 V27T lid pi Dt a 



Av E,, 



3 Av Yd 



(33) 



where we have used that m = (47r/3)pi nt a 3 and A = na 1 . What 
we do here is to compare f grow with the timescale for the radial 
drift given by 



fdrift = 



d\n r 



dt 



(34) 



Now we focus on the stage at which the radial drift velocity 
reaches the maximum value, i.e., Q.t s = 1. At this stage, the 
dust scale height is given by hd ~ {2a D /3) l ^ 2 h g according 
to Equation ^3), In addition, we set Ay » -\faoc s since the 
collision velocity between Q.t s = 1 particles is dominated by 
the turbulence-driven velocity. Using these relations and h g = 
cJQ., we can rewrite Equation ( 1331 as 



f gr°wlf!f,= l — — 




(35) 



where tK = 2n/Q. is the Keplerian orbital period. Thus, 
the growth timescale is shorter when the mass-to-area ratio 
m/A oc p int a is smaller. Note that f glow |nr,=i is independent of 
a D since both hd and Ay scale with yfao. By contrast, the 
drift timescale for Q.t s - 1 particles is 



fdriftlof,=l 



1 



40 



n 



4x 10" 



tK- 



(36) 



The ratio of the two timescales is 




(37) 



The ratio (f gr o W /fdrift)n/,=i determines the fate of dust growth 
at £lt s = 1. If (f g rowAdrift)nr,,=i is very small, dust particles 



grow beyond Qt s = 1 without experiencing significant radial 
drift; otherwise, dust particles drift inward before they grow. 
We expect growth without significant drift to occur if 



? drift 



n;,=i 



1 

30' 



(38) 



where the threshold value 1/30 takes into account the fact that 
fgrow is the timescale for mass doubling while the particles ex- 
perience the fastest radial drift over decades in mass. Below, 
we examine in what condition this requirement is satisfied. 

The ratio (fimtCt)nt,=i/^g depends on the drag regime at 
Qf s = 1. We consider the Epstein regime first. Using 
p g = Z 9 /( y/lnhg) and h g = c f /Q, one can rewrite Epstein's 
law as Qf ? = (7r/2)pj nt a/£ ff . Thus, for Qt s = 1, we have a 
surprisingly simple relation 



(Pmta)(it s =i _ 2 
Inserting this relation into Equation 



fgrow \ot,=l ~ 30 



0.01 



we obtain 



Ik- 



(39) 



(40) 



Hence, the growth condition (Equation (l38l ) is not satisfied 
for the standard disk parameters tj as 10~ 3 and Hd/^g 
0,0 1 . in am cement with the results of our and previous 
dBrauer et al.l 12008a!) simulations. Note that the right-hand 
side of Equation (f40b is independent of p; nt . Thus, the poros- 
ity of aggregates has no effect on the radial drift barrier within 
the Epstein regime. 

The situation differs in the Stokes drag regime. A similar 
calculation as above leads to 



and 



(Pinta)n/ I =i 



fgrow \at s =] ~ 60 



9 A 



mfp 



2n a\ 



0.01 



m s =i 



/Irnfp 

a\nt.=i 



(41) 



(42) 



Note that the growth timescale is inversely proportional to 
the aggregate radius, in contrast to that in the Epstein regime 
(Equation d40"l i) being independent of aggregate properties. 
Substituting Equations d36t and d42l into the growth condition 
(Equation (|38l), we find that aggregates break through the ra- 
dial drift barrier in the "deep" Stokes regime, a|nr,=iMmfp Z 
45. Unlike Equation d40l ). Equation ( l42l implicitly depends 
on pi n t through aaf s =i/^mfp (see below), so the porosity of 
aggregates does affect the growth timescale in this case. It 
is interesting to note that the speedup of dust growth occurs 
even though the maximum collision velocity is the same. In- 
deed, the collision velocity depends only on Qf s and is thus 
independent of the drag regime. We remark that Stokes' 
law breaks down when a becomes so large that the particle 
Reynolds number bec ome s much larger than unity, as alre ady 
mentioned in Section 12.21 We will show in Section 15.11 that 
this fact sets the minimum value of t smw \ci ts =\ to » 0.3f^; see 
Equation d47l ). 

The internal density of aggregates controls the growth 
timescale through the aggregate size a at Q.t s = 1 . For given 
Pint, one can analytically calculate a\nt„=i from Equations ( [391 
and fiTT ). Explicitly, 



a\nt s =i = 



2Zg_ 

TiPint 



(43) 
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Figure 11. Size a (upper panel) and growth timescale f al - ow (lower panel) 
of dust aggregates at Clt s = 1 as a function of orbital radius r for inter- 
nal densities pj nt = 1.4 gem" 3 (solid line), 10~ 2 g cirT 3 (dashed line), and 
10~ 5 g cirT 3 (dotted line). The MMSN with the dimensionless diffusion co- 
efficient ao = 10~ 3 is assumed for the disk model. The thick line in the upper 
panel indicates a = 9/i m f p /4, at which the drag law changes from the Epstein 
regime to the Stokes regime. The thick line in the lower panel shows the 
drift timescale tfom at Slt s = 1 (independent of p mt ). For p mt = 10~ 5 gem -3 , 
'eiow|oj s =i satisfies the growth criterion (Equation J381 ) at r < 1 AU . In 
reality, fg ro wln<,=l does not fall below the value given by Equation 1471 (thin 
dotted line) because of th e effect of the gas drag at high particle Reynolds 
numbers (see Section[5~lJ. However, this does not change the location where 
the growth condition is satisfied. 

for the Epstein regime, and 

3 / m„h a 
« ».=i = 75-TT74 A — 2 " 2 " (44) 

for the Stokes regime, where we have used A m f p = 

m g l(p g cr mo \) and p g - Z g /( y2nh g ). For fixed p int , a\ot s =i 
decreases with increasing r in the Epstein regime, but in- 
creases in the Stokes regime. The upper panel of Figure QT| 
plots a|n fj =i for three different values of the aggregate inter- 
nal density p mt . If dust particles grew into compact spheres 
(pint ~ 1 g cm -3 ), Epstein's law governs the motion of Q.t s = 1 
particles in almost entire parts of the snow region (r > 3AU). 
However, if dust particles grow into highly porous aggregates 
with pi nt ~ 10~ 5 g cirT 3 , the particles growing at r < 60 AU 
enter the Stokes regime before Qf s reaches unity. The lower 
panel of Figure [TT] shows the two timescales f gl0 wlnr,=i and 
fdriftlnr,=i as calculated from Equations d35l l and ( l36*l l. respec- 



Z g =10xMMSN, Z rf =0.01Z rf 



io 6 F 




10 1 io 2 
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Figure 12. Same as Figure ITT] but fo r a d isk 10 times heavier than the 
MMSN. The growth criterion (Equation J38I ) is satisfied at r < 25 AU for 
Pint = 10~ 5 g cm -3 and at r < 7 AU for pj nt = 10~ 2 g cirT 3 . 

tively. We see that compact particles with pj nt ~ 1 g cm~ 3 do 
not satisfy the growth condition (Equation (l38l l) outside the 
snow line, while porous aggregates with p; nt ~ 10 -5 g cirT 3 
do in the region r < 10 AU. These explain our simulation 
results presented in Section [3] 

Finally, we remark that a high disk mass (i.e., a high S„ with 
fixed Hd/^g) favors the breakthrough of the radial drift barrier. 
Figure[T2]shows the size a and the timescales f glow and fdrift at 
Qf. s = 1 for a disk 10 times heavier than the MMSN. We see 
that the growth condition (Equation ( TJST l) is now satisfied at 
r < 25 AU for p; nt = 10~" 5 g cm -3 and at r < 7 AU even 
for pi nt = 10~ 2 g cirT 3 . This is because a higher S 9 leads 
to a shorter A m f p and hence allows aggregates to reach the 
Stokes regime a/ Amfp » 1 at larger r or with higher p; nt (note 
that enhancement of Z 9 by a constant remains r\ and hence 
fdiiftlnr,=i unchanged). Interestingly, our porosity model pre- 
dicts that pintlnr s =i is independent of Y, g . In fact, substituting 
Equation (PHI with (Ai>)n, i= i » y[ar>c s and oc pi nt a 3 into 
Equation ( |28T i. we obtain the equation for pintlnf s =t that does 
not involve 1, g . 

5. DISCUSSION 

So far we have shown that the evolution of dust into highly 
porous aggregates is a key to overcome the radial drift bar- 
rier. On the other hand, in order to clarify the role of porosity 
evolution, we have ignored many other effects relevant to dust 
growth in protoplanetary disks. In this section, we discuss 
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how the ignored effects would affect dust evolution. 

5.1. Effect of the Friction Law at High Particle Reynolds 
Numbers 

In this study, we have assumed that the stopping time t s 
obeys Stokes' law whenever a > A m f v . In reality, Stokes' 
law applies only when the particle Reynolds number (the 
Reynolds number of the gas flow around the particle) Re p = 
2a\vd-v g \/v mo \ is less than unity, where \va-v g \ is the gas-dust 
relative velocity. When Re p > 1, i.e., the particle becomes so 
large and/or the gas-dust relative velocity becomes so high, 
the stopping time becomes depende nt on the particle veloc- 
ity (see, e.g., IWeidenschilling|fl977l) . In this subsection, we 
discuss how this effect affects our conclusion. 

In general, the stopping time at a » A m f p can be written as 



2m 



C D p g \vd - v g \A 



(45) 



where Cd is a dimensionless coefficient that depends on Re ; ,. 
Stokes' law, which applies when Re ; , «: 1, is given by Co = 
24/Re p . In the opposite limit, Re ; , » 1, the drag coefficient 
Co approaches a constant value (typically of order unity; e.g., 
Co ~ 0.5 for a sphere with 10 3 < Re ; , < 10 5 ), which is known 
as Newton's friction law. Thus, in the Newton regime, the 
stopping time depends on the particle velocity unlike in the 
Stokes regime. In this case, one has to calculate the stopping 
time and particle velocity simultaneously since the particle 
velocity in turn depends on the stopping time. 

In the previous sections, we have ignored the Newton 
regime to avoid the above-mentioned complexity. However, it 
is easy to calculate the growth timescale in the Newton regime 
for given Qf s , for which the gas-dust relative velocity can be 
known in advance. Below, we show that the Newton drag 
sets the minimum value of f grow |nf J =i (Equation (l35l l) for given 
orbital radius and internal density, which was not taken into 
account in Section [4] At the midplane, Equation (|45]> can be 
rewritten as Qf, = (2 y2~nlCo){Csl\Vd - v g \)m/(Y, g A), where 

we have used that p CJ = E g Q/( ^J2nc s ). When Qf, = 1, the 
gas-dust relative velocity is dominated by the dust radial ve- 
locity v r = -t]Vk, so we can set \v c { - v g \ « it\vk- Thus, at the 
midplane, we obtain a relation 



(pintfl)n/,=i ^ 3C, 
s 9 ~ 8 V2^ 



2n c s \0.5) c s 



where we have used that m/A - 4pi nt a/3. If Co reaches a con- 
stant, (pintfl)nf,=i/^9 no longer depends on aggregate proper- 
ties. Putting this equation into Equation ( l35l ). we have 



fgrowlfi/,= l ~ 0.3 



Zd/Zg\ l lC D \l 7]V K /C S 



0.01 



0.5 A 0.08 



Ik- 



(47) 



When Co = 24/Re p , Equation ( l47T i reduces to the equa- 
tion for the Stokes drag (Equation (l42l). where f gr0 wlnr,=i de- 
creases with increasing aggregate size a. However, when Re ; , 
becomes so large that Cd reaches a constant value, f gr0 wln( s =i 
no longer decreases with increasing a. Thus, we find that 
the Newton drag sets the minimum value of fgrow[nr s =i- For 
our disk model, in which ZdfZg = 0.01 and t]Uk/c s = 
0.08(r/5 AU) 1 / 4 , the minimum growth timescale is » 0.2- 
03(C D /0.5)t K at r * 3-10 AU. 

Since the Newton drag regime was ignored in our model, 
the growth rate of aggregates was overestimated there at high 



Re ; ,. As seen in the lower panel of Figure QT| the growth 
timescale f gro wln/ s =i for the p m - 10~ 5 g cirT 3 aggregates falls 
below the minimum possible value given by Equation Wf\ at 
r < 7 AU. This implies that dust growth is somewhat artifi- 
cially accelerated in our simulation presented in Section 13.21 
However, this artifact is not the reason why porous aggregates 
grow across the radial drift barrier in the simulation. Indeed, 
the drift timescale f grow |nr,=i is ~ 40/# at these orbital radii, 
and hence the minimum growth timescale still satisfies the 
condition for breaking through the drift barrier, Equation ( 1381 
(see Section |4). Thus, highly porous aggregates are still able 
to break through the radial drift barrier even if Newton's law 
at high particle Reynolds numbers is taken into account. 

In summary, we have shown that Newton's friction law 
(Cd ~ constant) at high particle Reynolds numbers sets a 
floor value for the grow timescale at Q f, = 1. In the nu- 
merical simulation presented in Section 13.21 the neglect of 
the Newton drag regime causes artificial acceleration of the 
growth of Qfj > 1 aggregates. However, comparison with the 
drift timescale shows that the floor value of f gro wlnf,=i is suffi- 
ciently small for dust to grow across Of , = 1 . Therefore, the 
deviation from Stokes' law at high particle Reynolds numbers 
has little effect on the successful breakthrough of the radial 
drift barrier observed in our simulation. 

5.2. Effects of Frictional Backreaction 

So far we have neglected the frictional backreaction from 
dust to gas when determining the velocities of dust aggregates 
(Equations ® and (TT8l). Here, we discuss the validity of this 
assumption. 

5.2. 1 . Effect on the Equilibrium Drift Velocity 

Frictional backreaction generally modifies the equilibrium 
velocities of both gas and dust. The equilibriu m velocities in 
the pre sence of the backreaction are derived by Ta naka et al.l 
(120051) for arbitrary dust size distribution. The result shows 
that the radial and azimuthal velocities v r and u^, - u'^ + vk of 
dust particles with stopping time t s are given by 



1 



1 + (Qf s ) 2 9 
Q^ 

2[l + (Qf s ) 2 ]' 



2Qf t 



1 + (Qf s ) 2 9A 
1 



1 + (Qf s ) 2 



where 



2T 



vg.r 



"94 



(l+X) 2 + Y 

1 +X 
"(1 +X) 2 + Y 2 



-T]V K , 



(48) 



(49) 



(50) 



(51) 



are the radial and azimuthal components of the gas velocity 
relative to the local circular Keplerian motion, respectively, 
and 

" Pd(m) T dm, (52) 



X 



I 

! 



p g 1 + (Qf s (m)) 2 

p d (m) Qf^(m) 
p g 1 + (Qf. s (m)) 2 



dm, 



(53) 



with pd(m) being the spatial mass density of dust particles per 
unit aggregate mass0 In the limit of X, Y — > 0, the gas ve- 

5 Equations j48M53l are equivalent to the "multi-species NSH solution" 
of Bai & Stone (2010a, their Equations (A4) and (A5)). 
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Figure 13. Radial and azimuthal velocities of gas (upper panel) and radial 
velocity of dust (lower panel) at r = 5 AU as a function of the weighted av- 
eraged mass (m) m . The solid black and gray curves in the upper panel show 
Vg tr and vF . = v gj f, - vk, respectively, obtained from the simulation includ- 
ing porosity evolution of aggregates and fractional backreaction from dust to 
gas. The dotted curves are the velocities when the fractional backreaction is 
neglected. 

locities approach v g<r — > and v' gi p — > -tjvk, and hence Equa- 
tions d48l i and (l49t reduce to Equations (O and (TTSt . respec- 
tively. Thus, the dimensionless quantities X and Y measure 
the significance of the frictional backreaction. As found from 
the integrands in Equations © and ( fT~8T >, the backreaction is 
nonnegligible when the local dust-to-gas mass ratio exceeds 
unity and the aggregates dominating the dust mass tightly cou- 
ple to the gas. 

To test the effect of fractional backreaction, we have also 
simulated porous aggregation using Equations (l48l and (|49l 
instead of Equations © and ( TT~8T > for the aggregate veloci- 
ties. However, it is found that the effect of backreaction is so 
small that the resulting dust evolution is hardly distinguish- 
able from that presented in Section [3] The upper panel of 
Figure [13] shows the temporal evolution of the gas veloci- 
ties v a r and v' , observed in this simulation as a function of 
the weighted average mass {m) m . We see that the observed 
gas velocities deviate at most only by 9 m s _1 » 0.17?/y^ 
from the velocities when the backreaction is absent (dotted 
lines). As a result of this, the inward velocity -v r of aggre- 
gates with m = {m) m is decreased only by 15 % even when 
Q.t s {(m) m ) « 1 (see the black solid curve in the lower panel 
of Figure [T3]l. The above result can be understood in the fol- 
lowing way. As found from the definitions of X and Y (Equa- 
tions (p2l and (1531). the effect of the backreaction is signifi- 
cant only when the density of dust coupled to the gas (Qf , < 
1) is comparable to or higher than the gas density. When 
Qr J ((m),„) < 1, the density of the coupled dust at the mid- 
plane is < I.d/h d \at,=i ~ Zd/(hg V»d) ~ (0.01/ VoD)p ;i ,mid ~ 
0.3p 9 , m id < p g ,mid, where p 9 , m jd is the midplane gas density and 



we have used that hd\m,=i ~ yfaBhg ~ 0.03/z 9 (Equation (O) 
and » 0.01 (the latter is true as long as Qt s ((m) m ) < 1). 

When Qt s ((m) m ) > 1, the dust density does exceed the gas 
density at the midplane, but the most part of the dust mass is 
now carried by decoupled (Q? s > 1) aggregates, which do not 
affect the gas motion^ Thus, the density of coupled dust is 
always lower than the gas density, and hence the backreaction 
effect is insignificant at all times. 

Furthermore, the effect on the differential drift velocity is 
even less significant, because the decreases in the inward ve- 
locities nearly cancel out. As an example, the gray solid and 
dotted curves in the lower panel of Figure [13] show the dif- 
ferential radial velocity Ay,, between aggregates of stopping 
times t s = t s ((m) m ) and Q.3t s ((m) m ) obtained from the sim- 
ulations with and without the backreaction, respectively. We 
see that the maximum values of \Av r \, which are reached when 
Q? s ((m) m ) « 0.7, differ only by 5%. Therefore, the frictional 
backreaction from dust to gas hardly affects the drift-induced 
collision velocity between dust aggregates. 

5.2.2. Streaming Instability 

The backreaction of dust on gas causes another 
phenomenon, the so- called streaming instability 
dYoudin & Goodman] 120051) . This means that the equi- 
librium gas-dust motion as described by Equations ©-(US) 
is unstable against perturbation. One important consequence 
of this instability is rapid clumpin g of marginally decoupled 
(Qt s ~ 1) dust particles (e.g.. pohansen & Youdinl 120071: 
Uohansen et alj|2007l; iBai & Stonel2010a) . The clumping pro- 
ceeds in a runaway manner (i.e., turbulent diffusion no longer 
limits the clumping) once th e dust density exceeds the gas 
density at the midplane (e.g., iJohansen & Youdinl 120071: se e 
also the analytic explanation of this by IJohansen et al.H2009l) . 
The runaway clumps could be eventually gravitationally 
bound and form 100 km sized planetesimals (IJohansen et al.l 
12007b . For more tightly coupled (Of 5 <K 1) particles, 
however, the clumping occurs only moderately unless the 
dust-to-gas s urface density ratio is high and/or the radia l drift 
speed is low (IJohansen et al.ll2009t IBai & Stonell2010bl) . This 
is also true for loosely coupled particles (Of, » 1) for which 
the interaction with the gas is weak. 

As seen in Section [3721 porous aggregates are able to reach 
£2fj ~ 1 in inner regions of disks. These aggregates likely trig- 
ger the streaming instability and can even experience runaway 
collapse. However, it is not obvious whether the clumps really 
experience the runaway collapse, since the growth timescale 
of the Clt s ~ 1 aggre gates can be as short as one orbital pe- 
riod (see Section lSTl i, which is compara ble to the growth time 
of the streaming instability at Q.t s = 1 dYoudin & Goodman! 
12005b . If the aggregates cross Qf s ~ 1 faster than the clumps 
develop, planetesimal formation will occur via direct colli- 
sional growth rather than gravitational instability. In order to 
address this issue, we will need to simulate coagulation and 
streaming instability simultaneously. 

5.3. Fragmentation Barrier 

In this study, we have assumed that all aggregate collisions 
lead to sticking. This assumption breaks down if the colli- 
sional velocity is so high that the collision involves fragmen- 
tation and erosion. If the mass loss due to fragmentation and 

6 Indeed, X and ¥ are insensitive to fif, 3> 1 particles because the factors 
l/[l+(fi(j) 2 ] » r s 2 andfi/,/[l + (fi/,) 2 ] ss r s l decrease faster than the spatial 
dust density pj oc Xj/hj oc f'' 2 increases (see Equation (3)) . 
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erosion is significant, it acts as an obstacle to pla netesimal for- 
mation (the so-called fragmentation barrier; e.g. jBrauer et alj 
l2008al) . Here, we discuss the validity and possible limitations 
of this assumption. 

Recent A^-body simulations predict that very fluffy aggre- 
gates made of 0.1 /mi sized icy particles experience catas- 
trophic disruption at collision velocities Ai> > 35 m s _1 
dWada et al .1120091) . If a large aggregate grows mainly through 
collisions with similar-sized ones (which is true in our simu- 
lations; see Figure|9]l, the collision velocity at Qt s as 1 is dom- 
inated b y the t urbulence-driven velocity Av, as 5v g as y[a~£,c s 
(Section [2.3.21 ). If the disk is optically thin and moderately 
turbulent {ao = 10~ 3 ) as in our model, the collision veloc- 
ity is as 21 m s _1 at r = 5 AU, so catastrophic disruption is 
likely insignificant for such collisions. However, if turbulence 
is as strong as ao = 10 -2 , the collision velocity at r — 5 AU 
and £lt s = 1 goes up to 67 m s . In protoplanetary disks, 
strong turbulence with ao > 10 2 can be driven by magne - 
torotational instability (MRI; e.g., iBalbus & Hawlevl[l998l) . 
If such strong turbulence exists, fragmentation becomes no 
more negligible even for icy aggregates. Besides, the colli- 
sion velocity can become higher than the above estimate when 
a large aggregate collides with much smaller ones, since the 
collision velocity is then dominated by the radial drift mo- 
tion. For example, the differential radial drift velocity be- 
tween an Qf j = 1 aggregate and a much smaller one is as high 
as ai tjuk ~ 56 m s _1 in optically thin disks. At such a high 
velocity, erosion by small aggregates can also slow down the 
growth of Q.t s ai 1 aggregates, although net growth might be 
possible (see, e.g.. lTeiser & Wurmll2009HTeiser et al.ll201 lh . 

On the other hand, resupply of small dust particles by frag- 
mentation/erosion has positive effects on dust growth. First, 
small dust particles stabilize MRI-driven turbulence because 
they efficiently capture ionized gas particles and thereby re- 
duce the electric conductivity of the gas (e.g., ISano et all 
I2000T) . This process generally leads to the reduction of the 
gas random velocity (and hence the reduction of turbulence- 
induced collision velocity), espe cially when the magneti c 
fields threading the disk are weak (lOkuzumi & Hirosell201 ll) . 
In addition, small fragments enhance the optical thickness of 
the disk, and thus reduce the temperature of the gas in the 
interior of the disk (given that turbulence is stabilized there). 
Since the radial drift velocity is proportional to the gas tem- 
perature, this leads to the reduction of the drift-induced col- 
lision velocity. In the limit of large optical depths, the gas 
temperatu re is reduced by a fa ctor us (h/r) 1 ^ 4 as 0.5 near the 
midplane dKusaka et al. 1970), resulting in the reduction of 
the drift-induced collision velocity to 28 m s _1 . These effects 
may help the growth of large aggregates beyond the fragmen- 
tation barrier. 

The size of monomers is another key factor. Al- 
though we have assumed monodisperse monomers of ao = 
0.1//m, the size of interstellar dust particles ranges from 
nanometers to microns. It is suggested both theoretically 
dChokshi et aT1ll993t iDominik & Tielenslll997l) and experi- 
mentallv dBlum & Wurm 1 120081) that the threshold velocity 
for sticking is roughly inversely proportional to ao. Thus, 
inclusion of larger monomers generally leads to the de- 
crease in the sticking efficiency. However, it is not obvious 
whether aggregates composed of multi-sized interstellar par- 
ticles are mechanically weaker or stronger than aggregates 
considered in this study. For example, if the monomer size 
distribution dno/dao obeys that of interstellar dust particles, 



OHo/aTogao °c a^ 2 dMathis et al.lll977T) . the total mass of 
the aggregates is dominated by the largest ones (mo °c a^ 

1 II 

and hence modno/dlogao oc a Q ). Nevertheless, the ex- 
istence of smaller monomers can still be important, since 
the binding energy per contact Ehreak is propo rtional to aJ 3 
dChokshi et all 19931 IDominik & Tielensll 19971) and hence the 
total binding energy tends to be dominated by the smallest 
ones (Zibreakflno/^log ao °c a^l 6 ). The net effect of multi- 
sized monomers needs to be clarified by future numerical as 
well as laboratory experiments. 

Another issue about the growth efficiency of icy aggre- 
gates arises from sintering. Sintering is redistribution of ice 
molecules on solid surfaces due to vapor transport and other 
effects. In this process, ice molecules tend to fill dipped sur- 
faces (i.e., surfaces with negative surface curvature) since the 
equilibrium vapor pressure decreases with decreasing the sur- 
face curvature. In an aggregate composed of equal-sized icy 
monomers , this process le ads to growth of the monomer con- 
tact areas dSironoll201 lH) and consequently to enhancement 
of the aggregate's mechanical strength such as F m \\. Sig- 
nificant growth of the contact areas could cause the reduc- 
tion of the aggregate's sticking efficiency since the dissipa- 
tion of the collision energy throug h internal rol ling/sliding 
motion could then be suppressed dSironol fr999). Further- 
more, if the monomers have different sizes, sintering leads 
to evaporation of smaller monomers (having higher positive 
curva ture), which m ay result in the breakup of the aggre- 
gate dSironol [201 lal) . Therefore, sintering can prevent the 
growth of icy ag gregates near th e snow line where sintering 
proceeds rapidlv. lSironol (1201 lbl) shows that the timescale of 
H2O sintering falls below 10 3 yr in the region between the 
snow line (3 AU) and 7 AU for the radial temperature adopted 
in our study. This is comparable to the timescale on which 
submicron-sized icy particles grow into macroscopic objects 
in this region (see Figure |7). However , if the disk is passive 
and optically thick dKusaka et alj|1970t) . no icy materials (in- 
cluding H2O an d CO2) undergo rapid sintering at r > 4 AU 
(Sirono 2011b). Moreover, the required high optical depth 
can be provided by tiny fragments that would result from the 
sintering-induced fragmentation itself. Consistent treatment 
of the two competing effects is necessary to precisely know 
the location where sintering is really problematic. 

To summarize, whether icy aggregates survive catastrophic 
fragmentation and erosion crucially depends on the environ- 
ment of protoplanetary disks as well as on the size distribution 
of the aggregates and constituent monomers. However, we 
emphasize that icy aggregates can survive within a realistic 
range of disk conditions as explained above. Indeed, the range 
is much wider than that for rocky aggregates, for which catas- 
trophic _chsnr£ti£n_o£cms^ 

m s" 1 dBlum & Wurm ||2508t IWada et al.ll2009t iGiittler et all 
120101) . In order to precisely predict in what conditions icy 
aggregates overcome the fragmentation barrier, we need to 
take into account the mass loss due to fragmentation/erosion 
and the reduction of collision velocities due to the resupply of 
small particles in a self-consistent way. This will be done in 
our future work. 

5 .4. Validity and Limitations of the Porosity Model 

Aggregates observed in our simulation have very low in- 
ternal densities. This is a direct consequence of the porosity 
model we adopted (Equation (fTBTl). Here, we discuss the va- 
lidity and limitations of our porosity model. 
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As mentioned in Section l2.3.11 our porosity change recipe 
at E[ mp > E m \\ is based on head-on collision experiments of 
similar-sized aggregates. In our simulation, dust growth is 
indeed domi nated by collision with similar-sized aggregates 
(see Section 13.2. It , so our result is unlikely affected by the 
limitation of the porosity model regarding the size ratio. By 
contrast, the neglect of offset collision may cause underesti- 
mation of the porosity increase, since the impact energy is 
spent for stretching rather than compaction a t offset collision 
(IWada et all [20071: iPaszun & Dominikl 120091) . If this is the 
case, then the breakthrough of the radial drift barrier can oc- 
cur even outside 10 AU. 

On the other hand, the formation of low-density dust aggre- 
gates is apparently inconsistent with the existence of massive 
and much less porous aggregates in our solar system. For 
example, comets, presumably the most primitive dust "aggre- 
gates" in the solar system, are exp ected to have mean interna l 
densities of p mt ~ 0.1 g cm 3 (e.g. JGreenberg & Hagelfl990l) . 
Since our porosity model does not explain the formation of 
such large and less porous "aggregates," there should exist 
any missing compaction mechanisms. 

One possibility is static compression due to gas drag and 
self-gravity. Although static compression is ignored in our 
porosity model, it can contribute to compaction of aggre- 
gates that are massive or decoupled from the gas motion. For 
relatively compact {p mt ~ 0.1 gem -3 ) dust cakes made of 
micron-sized Si02 particles, static com paction is observed 
to occur at static press ure > 100 Pa dBlum & Schraplerl 
l2004tlGtittler et all2009l) . By contrast, the static compression 
strength has not yet been measured so far for icy aggregates 
with very low internal densities (p ml «; 0.1 g cm). How- 
ever, for future reference, it will be useful to estimate here the 
static pressures due to gas drag and self-gravity. 

The ram pressure, the gas drag force per unit area, is given 
by -Pram = CoPg\vd — v g \ 2 /2, where Co is the drag c oeffic ient 
and \Vd~v g \ is the gas-dust relative speed (see Section lBTTl i. At 
Q.t s > 1, the gas-dust relative speed is approximately equal to 
tjvk- Thus, assuming New ton's drag law Co ~ 1 for Q.t s > 
1 aggregates (Section 15.1b , the ram pressure at Q.t s > 1 is 
estimated as 



Pram ~ Pg(rjV K ) ~ 10 



Pg 



10- 11 gem- 3 A 50 ms- 1 



Pa (54) 



independently of aggregate properties. Thus, if the static com- 
pression strength of our high porous aggregates is lower than 
10~ 5 Pa, compression of the aggregates will occur at Q.t s > 1 
due to ram pressure. 

The static pressure due to self-gravity is estimated from di- 
mensional analysis as 



Gm 2 



ur 



10 



10 



2/3 



Pint 



10- 5 g cm- 3 



4/3 



Pa. (55) 



For m ~ 10 10 g and p mt ~ 10 -5 g cuT 3 , which correspond 
to the Qi s = 1 aggregates observed in our simulation (Fig- 
ure [5), the gravitational pressure is much weaker than the 
ram pressure. However, since Pg ra v K m 2 ^, compression 
due to self-gravity becomes important for much heavier ag- 
gregates. For example, if pi nt ~ 10- 5 (m/10 I() gY 1 ^ g crrr 3 
as is for the Qt$ > 1 aggregates observed in our simulation, 
Pgrav exceeds P ra m at m ~ 10 g, which is comparable to the 
mass of comet Halley. Moreover, since P gr av K p 4/ ? ; gravi- 
tational compaction will proceed in a runaway manner unless 



the static compression strength increases more rapidly than 
Pgrav Thus, static compression due to self-gravity may be a 
key to fill the gap between our high porous aggregates and 
more compact planetesimal-mass bodies in the solar system. 

6. SUMMARY AND OUTLOOK 

We have investigated how the porosity evolution of dust 
aggregates affects their collisional growth and radial inward 
drift. We have applied a porosity model based on jV-body sim - 
ulations of aggregate collisions (ISuvama et alj|2008l 120121) . 
This porosity model allows us to study the porosity change 
upon collision for a wide range of impact energies. As a 
first step, we have neglected the mass loss due to collisional 
fragmentation and instead focused on dust evolution outside 
the snow line, where aggregates are mainly composed of ice 
and hence catastr ophic fragmentation may be insignificant 
(IWada et al. 2009]). Our findings are summarized as follows. 

1. Icy aggregates can become highly porous even if colli- 
sional compression is taken into account (Section [T2l . 
Our model calculation suggests that the internal density 
of icy aggregates at 5 AU falls off to 10~ 5 g cm- 3 by the 
end of the initial fractal growth stage and then is kept 
to this level until the aggregates decouple from the gas 
motion (Figure [TOt. Stretching of merged aggregates 
at offset collisions, which is not taken into account in 
our por osity model, could further decrease the inter nal 
density (IWada et al.ll2007HPaszun & Dominikll2009h . 

2. A high porosity triggers significant acceleration in col- 
lisional growth. This acceleration is a natural con- 
sequence of particles' aerodynamical property in the 
Stokes regime, i.e., at particle radii larger than the mean 
free path of the gas molecules (Section|4]i. The porosity 
(or internal density) of an aggregate determines whether 
the aggregate reaches the Stokes regime before the ra- 
dial drift stalls its growth. Compact aggregates tend 
to drift inward before experiencing the rapid growth, 
while highly porous aggregates are able to experience 
it over a wide range of the orbital radius (Figure [TT). 

3. The growth acceleration enables the aggregates to over- 
come the radial drift barrier in inner regions of the 
disks. Our model calculation shows that the break- 
through of the radial drift barrier can occur at orbital 
radii less than 10 AU in the MMSN (Figure |5}. A 
higher disk mass allows this to occur at larger orbital 
radii or higher internal densities (Figure [12}. The ra- 
dial drift barrier has been commonly thought to be one 
of the most serious obstacles against planetesimal for- 
mation. Our result suggests that, if the fragment ation 
of icy aggregates is truly insignificant (see Section l5T3l >. 
formation of icy planetesimals is possible via direct col- 
lisional growth of submicron-sized icy particles even 
without an enhancement of the initial dust-to-gas mass 
ratio. 

4. Further out in the disk, the growth of porous icy aggre- 
gates is still limited by the radial drift barrier, but their 
inward drift results in enhancement of the dust surface 
density in the inner region (Figure |6). This enhance- 
ment may help the core of giant planets to f orm within 
a disk lifetime (Kobayashi et al. 2 0101 12.0 fib . 
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We remark that the quick growth in the Stokes regime 
was also obs e rved in recent coagulation sim ulations by 
iBirnstiel et al.l (120101 see their Figure 1 1) an d IZsom et all 
(1201 11 see their Figure 3). IBirnstiel et al.l (120101) observed the 
breakthrough of the radial drift barrier only at small orbital 
radii (r < . 5 AU) since they assumed compact aggregation. 
IZsom et ail (1201 ll) found rapid growth of porous aggregates 
in the Stokes regime, but did not consider the loss of the dust 
surface density through radial drift. What we have clarified in 
this study is that porosity evolution indeed enables the break- 
through of the radial drift barrier at much lager orbital radii. 

The porosity evolution can even influence the evolution of 
solid bodies after planetesimal formation. It is commonly 
believed that the formation of protoplanets begins with the 
runaway growth of a small number of planetesimals due to 
gravitational focusing (e.g.. IWetherill & Steward [19 891) . The 
runaway growth requires a sufficiently high gravitational es- 
cape velocity u e sc = y/2Gm/a relative to the collision ve- 
locity. Since the escape velocity decreases with decreasing 
internal density (u esc oc m'^pV 6 ), it is possible that a high 
porosity delays the onset of the runaway and thereby affects 
its outcome. For example, a recent protoplan et growth model 
including co llisional fragmentation/erosion (Kobavas hi et alj 
120101 1201 11) suggests that planetesimals need to have grown 
to > 10 21 g before the runaway growth begins in order to en- 
able the formation of gas giant planets within the framework 
of the core accretion scenario (lMizunolll980l: iPollack et al.l 
[19961) . The size of the "initial" planetesimals can even 
deter mine the mass distribution of asteroids in the main 
belt dMorbidelli et alj|2009t lWeidenschilling|l201 11) . As we 
pointed out in Section 15.41 compaction of large and massive 
aggregates may occur through static compression due to gas 
drag or self gravity. To precisely determine when it occurs is 
beyond the scope of this work, but it will be thus important 
to understand later stages of planetary system formation. We 
will address this in future work. 
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